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Background: Many different shape degrees of freedom play crucial roles in determining the nuclear ground 
state and saddle point properties and the fission path. For the study of nuclear potential energy surfaces, it is 
desirable to have microscopic and self-consistent models in which all known important shape degrees of freedom 
are included. 


Purpose: By breaking both the axial and the spatial reflection symmetries simultaneously, we develop 
multidimensionally-constrained relativistic mean field (MDC-RMF) models. 


Methods: The nuclear shape is assumed to be invariant under the reversion of x and y axes, i.e., the intrinsic 
symmetry group is V4 and all shape degrees of freedom Bap with even u, such as $20, B22, 230, B32, Bao, .--, 
are included self-consistently. The single-particle wave functions are expanded in an axially deformed harmonic 
oscillator (ADHO) basis. The RMF functional can be one of the following four forms: the meson exchange or 
point-coupling nucleon interactions combined with the nonlinear or density-dependent couplings. The pairing 
effects are taken into account with the BCS approach. 


Results: The one-, two, and three-dimensional potential energy surfaces of 74°Pu are illustrated for numerical 
checks and for the study of the effect of the triaxiality on the fission barriers. Potential energy curves of even-even 
actinide nuclei around the first and second fission barriers are studied systematically. Besides the first ones, the 
second fission barriers in these nuclei are also lowered considerably by the triaxial deformation. This lowering 
effect is independent of the effective interactions used in the RMF functionals. Further discussions are made about 
different predictions on the effect of the triaxiality between the macroscopic-microscopic and MDC-RMF models, 
possible discontinuities on PES’s from self-consistent approaches, and the restoration of broken symmetries. 
Conclusions: MDC-RMF models give a reasonably good description of fission barriers of even-even actinide 
nuclei. It is important to include both the nonaxial and the reflection asymmetric shapes simultaneously for the 
study of potential energy surfaces and fission barriers of actinide nuclei and of those in unknown mass regions 
such as, e.g., superheavy nuclei. 


PACS numbers: 21.60.Jz, 24.75.+i, 25.85.-w, 27.90.+b 


I. INTRODUCTION years the Nilsson perturbed-spheroid parameter €2 


was usually adopted for numerical convenience [5]. 


The occurrence of spontaneous symmetry breaking 
leads to nuclear shapes with a variety of symmetries [1, 2]. 
A lot of nuclear phenomena are connected with the 
nuclear deformation, including small and large ampli- 
tude collective motions, e.g., the rotation and the fis- 
sion [3, 4]. The shape of a nucleus can be described by 
the parametrization of the nuclear surface or the nucleon 
density distribution. The multipole expansion of the nu- 
clear surface is usually used in mean field calculations, 


fore) À 
R(0,9) = Ro |1 + boo +X XO BYA (1) 
A=2 u=— À 


where (),,’s are deformation parameters. 

The majority of observed nuclear shapes is of the 
spheroidal form which can be described by the axial- 
quadrupole deformation parameter 629 though in early 
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The nonaxial-quadrupole (triaxial) deformation (22 (or, 
equivalently, y) manifests itself by the wobbling motion 
and chiral doublet bands [6-9], and it may also play im- 
portant roles in superheavy nuclei (SHN) [10]. 

The octupole shapes with À = 3 are predicted to ex- 
ist in nuclei in several mass regions [11]. The low-lying 
negative parity bands observed in actinides and some 
rare-earth nuclei are related to the reflection asymmetric 
(RA) shapes [12-17]. In addition, reflection asymmet- 
ric shell model calculations revealed that the observed 
low-energy 27 bands in N = 150 nuclei [18] are caused 
by the 832 deformation [19]. Indeed, strong Y32 cor- 
relations were found in some N = 150 isotones from 
multidimensionally-constrained covariant density func- 
tional theories (MDC-CDFT) [20]. 

Deformations of higher-order multipole with À > 3 are 
important to different extents. The hexadecapole defor- 
mation, 849 or €4, has been included in deformed mean 
field potentials since 1960s, see, e.g., Ref. [5]. The im- 
portant effects of the higher-order deformation {go or 
cs on the angular momentum alignments and dynamic 
moments of inertia in superheavy nuclei were also re- 
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vealed [21, 22]. 

The shape degrees of freedom are important not only 
for the ground states or small amplitude collective mo- 
tions, but also for large amplitude collective motions such 
as fission. Since the discovery of the nuclear fission, the 
description of the fission process has been a difficult and 
challenging task. The fission dynamics are mostly gov- 
erned by the barriers which prohibit the dissolving of 
the nucleus. In order to study the fission problem, one 
should have very accurate information about the fission 
barrier, i.e., the height and width, or, more precisely, the 
shape of the fission barrier [23, 24]. Particularly, to ex- 
plore the island of stability of SHN, it is more and more 
desirable to have accurate predictions of fission barriers 
of SHN. To now, many popular nuclear structure models 
have been employed to study nuclear fission barriers, in- 
cluding the macroscopic-microscopic (MM) models [25- 
29], the extended Thomas-Fermi plus Strutinsky integral 
(ETFSI) method [30], the Hartree-Fock or Hartree-Fock- 
Bogoliubov methods with the Skyrme force [31-38] and 
the Gogny force [39], and the CDFTs [31, 40-51]. 


Besides (29 which describes the elongation of a fis- 
sile nucleus and (649 which is relevant to the size of a 
neck, many other shape degrees of freedom are also cru- 
cial for determining the shape of fission barriers and the 
fission path. Let us take actinides as examples. Due 
to shell effects, actinide nuclei are characterized by a 
two-humped fission barrier [52]. It has long been known 
from the MM model calculations that the inner fission 
barrier is lowered by the nonaxial-quadrupole deforma- 
tion [53-55] and the outer one by the reflection asym- 
metric shape [56]. Later, the important roles played by 
the nonaxial-quadrupole deformation and the octupole 
deformation were confirmed in the nonrelativistic [57] 
and relativistic [58] density functional calculations, re- 
spectively. Therefore what is usually done is to consider 
the triaxial but reflection symmetric (RS) shapes for the 
inner barrier and axially symmetric (AS) but RA shapes 
for the outer one [39, 59, 60] though in several publica- 
tions, both the nonaxial and the octupole deformations 
are included [61]. 


In recent years, the nuclear CDFT has been very suc- 
cessful in the description of both ground states and ex- 
cited states of the nuclei ranging from light to superheavy 
regions. The CDFTs have also been used to study the 
PES’s and the fission barriers of heavy and superheavy 
nuclei [31, 40-51]. We have developed MDC-CDFTs by 
breaking both the axial and reflection symmetries simul- 
taneously [46-48]. In these MDC-CDFTs, all shape de- 
grees of freedom bap, with even u, e.g., B20, B22, B30, P32, 
Bao, ---, are included self-consistently. The covariant 
density functional can be one of the following four forms: 
the meson exchange or point-coupling nucleon interac- 
tions combined with the nonlinear or density-dependent 
couplings. For the particle-particle channel, either the 
BCS approach or the Bogoliubov transformation has 
been implemented. For convenience, we name the MDC- 
CDFT with the BCS approach for the pairing as the 


MDC-RMF models and those with the Bogoliubov trans- 
formation as the MDC-RHB models. Due to the heavy 
computational burden, MDC-RHB models have not been 
used to do multidimensionally constrained calculations 
yet (here “multi” means three or more). The MDC-RMF 
models have been used to study the PES’s of actinide nu- 
clei in the (820, 822, 830) deformation space and it was 
found that the triaxiality also plays an important role 
upon the second fission barriers [46]. In this paper, we 
will present the detailed formulae for MDC-RMF models 
and some results of actinide nuclei. 

The paper is organized as follows. The formalism of 
MDC-RMF models will be given in Sec. II. In Sec. III we 
present the numerical details and results on the PES’s of 
the actinide nuclei. A summary is given in Sec. IV. 


II. FORMALISM OF MDC-RMF MODELS 


In the RMF model, a nucleus is treated as a compos- 
ite of nucleons interacting through exchanges of mesons 
and photons [42, 62-68]. The effects of mesons are de- 
scribed either by mean fields or by point-like interactions 
between the nucleons [69, 70]. Meanwhile, the nonlinear 
coupling terms [71-73] or the density dependence of the 
coupling constants [74, 75] are introduced to give correct 
saturation properties of nuclear matter. Accordingly, the 
form of the covariant density functional can be one of the 
following four: the meson exchange or point-coupling nu- 
cleon interactions combined with the nonlinear or density 
dependent couplings. Most of the computational efforts 
are devoted to solving the Dirac equation and the calcu- 
lation of various densities; this is common for all these 
RMF models. In this section, we mainly present the 
formalism of the RMF model with the nonlinear point- 
couplings (NL-PC). 

The starting point of the RMF model with the NL-PC 
is the following Lagrangian, 


ee w(iy,0" Laer 


Myy Liin Lui Leoou;, (2) 


where 
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are the linear coupling, nonlinear coupling, derivative 
coupling, and the Coulomb part, respectively. M is the 
nucleon mass, «S, Av, ATS, ATV, ps, YS, VV, ds, ôv, 
org, and ôry are coupling constants for different chan- 
nels, and e is the electric charge. ps, Prs, py, and pry 
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are the isoscalar density, isovector density, time-like com- 
ponents of isoscalar current, and time-like components 
of isovector current, respectively. The densities and cur- 
rents are defined as 


ps=0b, prs = br, 
pv =Y, Prv = Y7. (4) 
Starting from the Lagrangian, using the Slater deter- 
minants as trial wave functions, and neglecting the Fock 
term as well as the contribution to the densities and cur- 
rents from the Dirac sea (the no sea approximation), we 


can derive the RMF equation with the variational prin- 
ciple, 


hbk (r) = exe (1), (5) 


where 


h=a-p+B|M+S(r)]+V(r), (6) 


is the single-particle Hamiltonian and 


S = asps t+arsprs-7+ Bspst+ysps 


+ dsAps + ôrsAPrs -7, (7) 
V = avpv +arverv T+ Wwpvypv 
+ dvApy + ôrvAfrv `T, (8) 


are the scalar and vector potentials, respectively. In the 
present work we suppose that the nuclei in question are 
invariant under the time-reversion, which means that all 
the time-odd or vector components of the currents and 
the potentials vanish. In this case the single-particle 
Hamiltonian has the time-reversal symmetry which sim- 
plifies the calculation. 

It is customary to solve the deformed RMF equations 
by expanding the auxiliary single-particle wave functions 
in a complete basis, e.g., the harmonic oscillator (HO) 
basis [76, 77] or the Woods-Saxon (WS) basis [78, 79]. 
By using a basis in a two-centre harmonic oscillator po- 
tential, a reflection asymmetric relativistic mean field 
(RAS-RMF) approach has been developed [80] and used 
to study PES’s of even-even !46-156Sm in which the 
important role of the octupole deformation on shape 
phase transitions was found [81]. However, in our MDC- 
CDFTs, the single-particle wave functions and various 
densities are expanded in an axially deformed harmonic 
oscillator (ADHO) basis. The ADHO basis consists of 
the eigenstates of the Schrödinger equation [76, 77, 82], 


h2 
—5q7 tP) Palro) = Ea®alro), (9) 
where r = (z, p) with p = yx? + y? and 
1 
Vg(z, p) = 5M (wie F wee"); (10) 


is the ADHO potential and w, and wy, are the oscillator 
frequencies along and perpendicular to the symmetry z 
axis, respectively. The solution of Eq. (9) reads 


Palro) = Cagn. (2) Ran (O) Xs (a), (11) 


where @n, (z) and Rẹ! (p) are the harmonic oscillator 
wave functions, 


1 1 z \ 2 
on. (2) = Je it fen (=) e z, (12) 
1 2n! 2 (5) ai 
RY (p) = =| | rll © ) e 08) 
“e (P) bp V (mp + |m])! (£ ne bs — 


Xs, is a two-component spinor and Cy is a complex 
number introduced for convenience. Harmonic oscilla- 
tor lengths b, and bp are related to the frequencies by 
bz = 1//Mw, and bp = 1/,/Mw,. The corresponding 
eigenenergy Ey = w,(2n, +|mi| +1) +wz(nz +1/2) and 
the major quantum number Na = 2n, + |mi| + nz. 
These basis states are also eigenstates of the z com- 
ponent of the angular momentum j, with eigenvalues 
Ka = m + ms. For any state ®a(ro), the time- 
reversal state is defined as ®a(rao) = T®,(ro), where 
T= ioy K is the time-reversal operator and K is the 
complex conjugation operator. Apparently, we have 
Ka = -Ka and ta = Ta, where 7g = +1 is the 
parity. The deformation of the basis (basis is defined 


through the relations w, = wo exp (-1 /5/4T basis) and 


Wp = wo exp (/5/16r bia) where wo = (wzw?) is 
the frequency of the corresponding spherical oscillator. 

These basis states form a complete set for expanding 
any two-component spinors. For a Dirac spinor with four 
components, 


neo- (Fofana) w 


where the sum runs over all the possible combination 
of the quantum numbers a = {n7, np; Mı, Ms} and fe 
and g® are the expansion coefficients. In practical cal- 
culations the summations in Eq. (14) have to be trun- 
cated. Following Ref. [83], for the large component of 
the Dirac wave function, the states satisfying [n./Q» + 
(2n, + |mi|)/Q,] < Ny are included in the expansion 
where Q, = max(1,b,/bo) and Q, = max(1,b,/bo) are 
constants related to the oscillator lengths bo = 1/V Mwo, 
bz, and bp. For the expansion of the small component, 
the truncation is made up to Ng = Ny + 1 major shells 
in order to avoid the spurious states [76]. 

If a nucleus is invariant under the rotation around the 
symmetry z axis and the spatial reflection, the angular 
momentum projection on the z-axis and the parity are 
conserved and the RMF equation (5) can be decomposed 
into blocks characterized by the quantum numbers Ka 
and Ta. Usually only half of the space with Ka > 0 is 
considered due to the time-reversal symmetry. 

Now let us turn to the general case that the axial sym- 
metry as well as the spatial reflection symmetry are bro- 
ken. Since we still expand the spinors in the ADHO 
basis, components with different K and 7 are mixed to- 
gether; thus, we must diagonalize a larger single-particle 
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Hamiltonian matrix with non-zero matrix elements be- 
tween two basis states with different K and m. Neverthe- 
less, even in this case, we still have one symmetry oper- 
ator that makes the Hamiltonian matrix block-diagonal. 
Due to the axial symmetry of the basis, it is convenient 
to introduce the simplex operator Ê = ie Note 
that for a fermionic system with a half-integer spin, S 
is a Hermitian operator and S$? = 1. This operator cor- 
responds to the rotation by m around the z axis, thus 
leaving the nuclear mean field invariant. The eigenvalue 
of S, S, is also a good quantum number for the basis, 
S®, = SPa = (—1)e—!/26,, which means that the ba- 
sis ®, with Ky = +1/2, —3/2, +5/2, —7/2, ... span 
the subspace with S = 1, while their time-reversal states 
span the one with S = —1. Note that now the blocks with 
K = +1/2, —3/2, +5/2, —7/2,... are mixed. Remem- 
ber that for a nucleus with the time-reversal symmetry, 
only the basis with S = 1 are used in the expansions; for 
such a basis state, we set Ca = 1 [cf. Eq. (11)] when 
expanding the large component and Ca = i for the small 
one. The basis states with S, = —1 are obtained by 
simply applying 7 on those with Sa = 1. Furthermore, 
for systems with the time-reversal symmetry, it is only 
necessary to diagonalize the matrix with S = 1 and the 
other half is obtained by a time-reversal operation on the 
obtained single-particle wave functions. 

For deformed nuclei with the V4 symmetry, we expand 
the potentials V(r), S(r), and the densities in Eq. (4) in 
terms of the Fourier series, 


—inje 


f(p, 9, 2 


=> fule,z = eplin). (15) 


p=— co 


Applying the symmetry conditions, it is easy to see that 
fu = fă = fa and fn = 0 for odd n. Thus the expansion 
(15) can be simplified as 


f(e, 9,2) = Fz olny), 


fol P; Z 


Fath P, Z 


(16) 
where 
27 


fo(p,2) _ dy f(p, p, 2) (17) 


Vir 


1 27 
,z) = — d , 9, z) cos(2nyp), 18 
fulor2)=—e | detlo z) cos(2no), (18) 
are real functions of p and z. The details for calculat- 
ing the matrix elements of the Dirac Hamiltonian and 
various densities and their derivatives are given in the 
Appendixes. 

For open shell nuclei the pairing interaction becomes 
crucial and must be included. It has been shown that fis- 
sion barriers depend very much on the form and strength 
of the effective pairing interactions [84]. Several methods 
have been developed to treat the pairing effects, e.g., the 
BCS approach, the Bogoliubov transformation, and the 
particle number conserving method [85-89]; all of them 


have been used in the study of PES’s and fission barri- 
ers. Since we use the BCS approach in our MDC-RMF 
calculations, we only show the formulas for the BCS ap- 
proximation. In the particle-particle channel, the gap 
equation reads [2], 


An= 05 r Aw (19) 


kkk’ k! = 7? 
e0? Ve + Ai 


where €, = €k — A and J is the Fermi energy. The total 
pairing energy is 


Ps aM i Perte. (20) 


In our models, either the 6-force or the finite-range sep- 
arable force [90] is implemented. 

To obtain a PES one can perform a constraint calcula- 
tion which is equivalent to adding an external potential 
during the iteration [2]. The quadratic constraint method 
is usually used, 


1 
E' = Erm +X 5 Oan (Qan — mu, (21) 
Au 


where C),, is the spring constant and m),,’s are desired 
moments. With this method the calculation always con- 
verges to a deformation point on the PES other than the 
desired one. To overcome this shortcoming and to get a 
PES with equally distributed points, we use a modified 
linear constraint method. The Routhian reads, 


1 
E’ = EgmF + 5 gO ru (22) 
Au 


where the variables C),,’s change their values during the 
iteration through the following relation, 


ole = E $ kau (e0 - Bus) (23) 


where (),, is the desired deformation, kau is a constant, 


and cm is the value at the nth step. This constraint 
method works well in our multidimensionally-constrained 
calculations. 

The total energy of the nucleus reads 


Etotal = per {Sieh (a “p+ BM) Vk 
k 


1 Dp 2 1 2 1 2 
+ 5asP5 + gov PV F JOTSPTS F a ATV PTV 


1 1 1 
| 3Bsps taia + FweV 


1 1 
+ zôsesAps + 5ovpvApy 


1 1 1 
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where Ee.m. is the center of mass correction. Depending 
on the effective interactions used in the RMF functional, 
Eom. can be calculated either in the oscillator approxi- 
mation, 


iws -; x 41A1/3 MeV, (25) 


or from the quasi-particle vacuum, 


(P?) 


Eom. = OMA’ 


(26) 
where P is the total linear momentum and A is the nu- 
clear mass number. 

The intrinsic multipole moments are calculated from 
the density by 


Qui = | Prove) Yal), (27) 


where Y), (Q) is the spherical harmonics. The deforma- 
tion parameter {, is obtained from the corresponding 
multipole moment by 


4T 
Ban = IN N (28) 


where R = 1.2 x A13 fm is the nuclear radius and N is 
the number of protons, neutrons, or nucleons. 


II. RESULTS AND DISCUSSIONS 
A. Numerical details 


In this section we present the numerical details and 
some illustrative calculations of MDC-RMF models. The 
potentials and densities are calculated in a spatial lattice 
[cf. Eqs. (15), (17), and (18) and Appendix B] in which 
mesh points in the p and z directions are designed in a 
way that the Gaussian quadrature can be made and those 
for the azimuthal angle ¢ are equally distributed. Since 
we keep the mirror reflection symmetry with respect to 
the x = 0 or y = 0 planes, only mesh points with positive 
x and y are considered. For the azimuthal angle ¢, more 
than 10 mesh points for light nuclei and about 20 mesh 
points for heavy nuclei are used, which are enough for 
most of practical applications. In two special cases the 
number of mesh points can be further reduced: (1) for 
axially symmetric nuclei the azimuthal degree of freedom 
vanishes and (2) for reflection symmetric nuclei the mesh 
points with z < 0 can be omitted. With these choices 
the values of the localized fields and potentials in the 
full lattice space can be simply obtained by symmetry 
transformations such as rotations or the spatial reflec- 
tion. The Coulomb field must be treated carefully due 
to its long-range nature. In the pairing channel, we use 
a density-independent 6 force with a smooth cutoff. The 
pairing strength parameters are Vp = —349.5 MeV fm? 
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FIG. 1. (Color online) The potential energy curve of 7*°Pu 
from MDC-RMF calculations with different truncations of the 
ADHO basis. Both triaxial (TA) and reflection asymmetric 
(RA) deformations are allowed. The results calculated with 
Ny = 16, 18, and 20 are depicted by dashed, dotted, and solid 
curves, respectively. The four sub-figures show the detailed 
structure of the potential energy curve near the ground state 
(A), the inner barrier (B), the isomeric state (C), and the 
outer barrier (D). The results calculated with Ny = 22 (dot- 
dashed curves) are also included in the sub-figures. The width 
of each sub-figure is 0.1 and the height is 1 MeV. 


and V, = —330.0 MeV fm? which are obtained by fitting 
the average pairing gaps [91, 92]. 

The calculated physical observables should converge 
as the truncation Ny — oo. In Fig. 1 we show the po- 
tential energy curve of 74°Pu calculated with different 
truncations, Nr = 16, 18, and 20. The effective interac- 
tion PC-PK1 [92] is used; accordingly, Ee.m. is calculated 
with Eq. (26). As is found in most of earlier calculations, 
the results show a typical two-humped structure. To see 
more clearly the truncation errors we have amplified the 
figure near four important points, i.e., the ground state, 
the top of the inner barrier, the isomeric state, and the 
top of the outer barrier. In these four sub-figures, the re- 
sults with Ny = 22 are also shown and we can investigate 
in more details the convergence properties of our model. 
When Np increases from 16 to 20, the binding energy 
changes differently for different points in the potential 
energy curve; the largest changes are near the ground 
state and about 0.3 MeV. Around the ground state and 
the second minimum, 74°Pu is axially deformed; the en- 
ergy obtained from calculations with Nr = 20 and 22 
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FIG. 2. (Color online) Mean field energy (i.e., Ec.m. is not in- 
cluded) of *4°Pu calculated by using basis with different basis 
deformation basis and different Ny. The results calculated 
with Ny = 16, 18, and 20 shells are denoted by red diamonds, 
green squares, and blue dots, respectively. The deformation 
is constrained to 620 = 0.3, 1.3, and 2.0, respectively. 


are almost the same. This means a good convergence; in 
the present work, we use Ny = 20 in axially symmetric 
calculations. Around the two fission barriers, the triaxial 
deformation is very important. One finds that for the 
inner barrier, the results from Ny = 20 and 22 are also 
almost identical; for the outer one, the difference between 
the barrier heights from calculations with Ny = 20 and 
22 is about several tens keV. It is very time consuming 
to make calculations with both axial and reflection sym- 
metries broken. In our systematic calculations presented 
in this paper, Ny = 20 is used around the ground state 
and the first saddle point, Nf = 16 is used around the 
fission isomer and the second saddle point. As a result, 
the inner barrier height is described with an accuracy 
of ~ 0.15 MeV and the outer one with an accuracy of 
~ 0.4 MeV. If in future calculations, Ny = 20 can also 
be used around the second barrier, the accuracy for its 
height should be within 0.2 MeV. 

Next we show how the results depend on the deforma- 
tion of the ADHO potential which is used to generate 
the ADHO basis; for brevity, we will call it the basis de- 
formation and label it with Spasis. In Fig. 2 we depict 
the calculated mean field energy of 74°Pu as a function 
of the basis deformation. Note that Ee.m. is not included 
in Fig. 2. In principle, if the basis space is complete, 
the results should not change when the basis deformation 
changes. Near the ground state, G29 = 0.3, the calculated 
energies are rather stable against the basis deformation. 
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FIG. 3. (Color online) Mean field energy (i.e., Ec.m. is not 
included) of *°Pu calculated with different truncations on 
large (upper panel) and small component (lower panel). In the 
calculations for the upper panel, Ng are fixed to 23, while the 
calculations for the lower panel are preformed with Ny = 10. 
The deformation is constrained to B29 = 0.3. 


Furthermore, the results with Ny, = 16, 18, and 20 al- 
most coincide with each other. This conclusion holds also 
for the second barrier, G29 = 1.3, except that two points 
with small basis deformations and Ny = 16 are very high. 
For even larger deformations, 329 = 2.0, the results with 
Ny = 16 deviate a bit from those with Np = 18 and 20; 
the differences between the results with Ns = 18 and 20 
are still very small. In the present work, the basis defor- 
mation is chosen in the following way: pasis = 620 for 
B20 < 0.3 and pasis = 820/2 for 0.3 < B20 < 2.0. 

When the basis is not complete, the calculated single- 
particle energies contain the contributions from the other 
levels in both the Fermi sea and the Dirac sea. This 
causes the “variational collapse” problem [93, 94]. As 
shown in Fig. 3, if we fix the truncation for the small com- 
ponent N, and increase Ny alone, the nucleus becomes 
more and more bound. Note that we have to set Nf < Ng 
to prohibit the occurrence of the spurious states, as we 
mentioned earlier. However, if we fix the truncation for 
the large component Ny and increase N, alone, the nu- 
cleus becomes less and less bound. Thus, as more basis 
states are included in the expansion, the binding energy 
may not change monotonically; this is different from the 
nonrelativistic calculations. 

In the calculations, the oscillator length or, equiva- 
lently, the frequency of the oscillator potential hwo for 
the ADHO basis should also be chosen carefully. The 
dependence of the binding energy of 74°Pu on hwo has 
been investigated in detail. As shown in Fig. 4, when 
R = fiwo/(41A~!/3 MeV) increases from 0.8 to 1.2, the 
binding energy of ?4°Pu in the whole potential energy 
curve varies by less than 1 MeV (0.05% of the absolute 
value). Moreover, when fiwo is around 41A~!/3 MeV, 
the binding energy changes very slightly. So, we set 
hwo = 41A7!/3 MeV in all calculations. 


201707.00926v1 


chinaXiv 


£ 240s [oct R=0.8 PC-PK1 
1800 Pu ee m 
a B —— R=1.0 
> --—-- R=1.1 a 
= -1805 Gdi 
| aro 
-1810 Fa 


- B: Inner 


barrier 


C: Fission 
isomer 


A: Ground Ag 
© state 7%" 


a 0.1 = B 


FIG. 4. (Color online) The potential energy curve of 7°Pu 
calculated by using basis with different oscillator frequency 
hwo. The axial symmetry is imposed and Ny = 20 is used. 
The results calculated with R = hwo/(41A~‘/? MeV) = 0.8, 
0.9, 1.0, 1.1, and 1.2 are depicted by dotted, dashed, solid, 
dot-dashed, and dash-dot-dotted curves, respectively. The 
four sub-figures show the detailed structure of the potential 
energy curve near the ground state (A), the inner barrier (B), 
the isomeric state (C), and the outer barrier (D). The width 
of each sub-figure is 0.1 and the height is 2 MeV. 


B. Three-dimensional PES of 74°Pu 


The double-humped fission barriers of actinide nuclei 
are usually used as a benchmark for theoretical models, 
see, e.g., Refs. [35, 40, 43, 44, 57, 58, 88, 95-101]. Asa 
nucleus evolves from the ground state to the fission con- 
figurations, many shape degrees of freedom play impor- 
tant and different roles in determining the heights of the 
inner and outer barriers. It has long been known that the 
inner barrier is lowered when the triaxial deformation is 
allowed [53-55], while for the outer barrier, the reflection 
asymmetric shape is favored [56]. We have shown that 
the reflection asymmetric outer barrier may be further 
lowered by including the nonaxial shape degrees of free- 
dom [46]. In this section we show some results of multi- 
dimensional PES’s for actinides. We will first present de- 
tailed results about the three-dimensional PES of 74°Pu, 
then display systematic results of the actinide nuclei. In 
these calculations, the effective interaction PC-PK1 [92] 
is used. 


In Figs. 5 and 6 is shown the calculated three- 


Í 


dimensional PES of 74°Pu. In these calculations, both 
non-axial and reflection asymmetric shapes are allowed 
and Ny = 16 is used. In each sub-figure we fix the value 
of 829 and display the energy as a function of 822 and B39. 
In other words, we have made three-dimensional con- 
straints on the corresponding multipole moments. Note 
that other shape degrees of freedom §), with even p, 
e.g., 832, Bao, B4s2, Baa, P50, ---, are also included in 
the calculations self-consistently. The MDC-RMF equa- 
tions are solved for each point on the deformation lattice 
(B20, B22, 330) in which B20 runs from 0.25 to 1.70 with 
a step size of 0.05, 622 from 0 to 0.25 with a step size of 
0.01, and 830 from 0 to 0.50 with a step size of 0.05. The 
points with 822(830) < 0 are obtained through the rela- 
tion E620, 822, 830) = E(b20, |822]; |830|). That is, for 
each sub-figure, 26 (for $22) x 11 (for 639) = 286 points 
are calculated, and there are totally 286 x 30 (for 620) 
= 8580 points in Figs. 5 and 6. This deformation lattice 
covers the shape space of the most interest for 74°Pu, 
from the ground state to the isomeric state and fission 
configurations. 


Now we examine the first three sub-figures with Boo = 
0.25, 0.30, and 0.35. It is clear that the ground state with 
B29 ~ 0.3 is both axially symmetric and reflection sym- 
metric, though it is a little soft against the octupole dis- 
tortion. When the nucleus is stretched by the quadrupole 
constraining potential, it becomes softer against the tri- 
axial as well as the octupole distortions. From the sub- 
figures with B29 = 0.40 ~ 0.65 one finds that two sym- 
metric minima with non-zero triaxial deformation 622 ap- 
pear and the corresponding fission paths are much more 
favored than the axially symmetric one, which is con- 
sistent with earlier calculations [45]. Although 639 = 0 
for all minima in these sub-figures, the softness against 
the octupole distortion changes as the nucleus is elon- 
gated. The inner fission barrier locates near the defor- 
mation 29 ~ 0.60. For the last six sub-figures in Fig. 5 
with 629 = 0.70 ~ 0.95, the situation becomes much sim- 
pler, where the nucleus becomes axially symmetric and 
reflection symmetric again. Nevertheless, the trend to 
become softer against the octupole distortion can be seen 
from these sub-figures, which indicates that the reflection 
asymmetric shape becomes more relevant. 


As the deformation $29 becomes larger, the second 
minimum and the second saddle point of the PES ap- 
pear. From the last three sub-figures of Fig. 5 and the 
first three sub-figures of Fig. 6 we see that, the nucleus 
keeps reflection symmetric near the fission isomeric state 
with Boo ~ 0.95, but it becomes softer against the 830 
distortion. Note that the scales of the $39 axes are dif- 
ferent in Figs. 5 and 6. At Gg9 = 1.15, two minima corre- 
sponding to reflection asymmetric shapes appear. Here 
the effect of the nonaxial deformation is not apparent, 
but along the 622 direction the PES becomes softer. At 
B20 = 1.2 the energy of the reflection asymmetric shape is 
lower by about 1 MeV than that of the reflection symmet- 
ric one. Interestingly, when $29 increases further, around 
the top of the second barrier, each reflection asymmet- 
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FIG. 5. (Color online) Sections of the three-dimensional potential energy surface, E = E (620, 822,830), of *4°Pu around the 
ground state, the inner barrier, and the fission isomer from MDC-RMF calculations. In each sub-figure the energy is shown as 
a function of the deformation parameters 622 and 639 when (20 is fixed at a certain value. The energies are normalized with 
respect to the binding energy of the ground state. The contour interval is 1 MeV. 


ric minimum splits into two minima with non-zero ĝə22. 
This is the lowering effects of triaxiality on the outer bar- 
rier found in Ref. [46]. Around the second barrier, the 
largest energy gain due to the triaxial distortion is about 
1 MeV. The nucleus becomes axially symmetric again 
when {29 > 1.6. 

From the above discussions, we can draw the following 
conclusions for 74°Pu: (1) Both the ground state and the 
fission isomeric state are axial and reflection symmetric; 
(2) Around the first fission barrier it assumes triaxial 
and reflection symmetric shapes; (3) Around the second 
fission barrier both triaxial and octupole deformations 
are important. 


C. PES’s of even-even actinide nuclei around two 
fission barriers 


The self-consistent three-dimensionally-constrained 
calculations are very time-consuming, we have only per- 
formed such calculations for 74°Pu. From this benchmark 


study we learned many experiences about the important 
roles played by various shape degrees of freedom in differ- 
ent regions of the deformation space, including the con- 
clusions listed in the end of Sec. II B. These experiences 
are used in a systematic study of even-even actinide nu- 
clei. 


Since around the inner barrier an actinide nucleus as- 
sumes triaxial but reflection symmetric shapes, we can 
make a one-dimensional constraint calculation with the 
triaxial deformation allowed and the reflection symmetry 
imposed. In Fig. 7 we present potential energy curves of 
even-even actinide nuclei around the ground state and the 
first fission barrier from MDC-RMF calculations. The 
axially symmetric results are displayed by dotted curves 
and those from triaxial calculations are shown by solid 
curves. The empirical values of fission barrier heights are 
taken from Ref. [102] and shown with red dots. We also 
list in Table I the heights of the first and second fission 
barriers of actinide nuclei from MDC-RMF calculations 
compared with the empirical values. It is clearly seen 
that the triaxial deformation lowers the inner barrier of 


FIG. 6. (Color online) Sections of the three-dimensional potential energy surface, E = E(620, 822,830), of °° Pu around the 
fission isomer and the outer barrier from MDC-RMF calculations. In each sub-figure the energy is shown as a function of the 
deformation parameters 622 and $309 when £20 is fixed at a certain value. The energies are normalized with respect to the 
binding energy of the ground state. The contour interval is 1 MeV. 


these actinide nuclei by about 1 ~ 4 MeV and this is con- 
sistent with Ref. [45]. The agreement of our calculation 
results with the empirical ones is good for most of the 
nuclei studied here with exceptions in the two thorium 
isotopes and ?38:240U, We have discussed possible reasons 
for these discrepancies in Ref. [46]. 

It is more complicated to determine the second fission 
barrier height because more shape degrees of freedom be- 
come important and there are often two or more fission 
paths. What we have done is the following [46]: (1) The 
axial symmetry is assumed and a two-dimensional con- 
straint calculation in the (820, 830) plane is made; (2) 
From the two-dimensional calculation the lowest fission 
path Blest (829) is approximately identified; (3) Along 
this fission path, a one-dimensional (29-constrained cal- 
culation is performed with the triaxial and octupole de- 
formations allowed and at each point with (29 fixed, 
the initial deformations are taken as Gi = 0 and 

ini. = Blowest (B59): (4) In this one-dimensional poten- 
tial energy curve, the second saddle point is located and 
the second barrier height is extracted. 


Potential energy curves of even-even actinide nuclei 
around the fission isomer and the second fission bar- 
rier from MDC-RMF calculations are shown in Fig. 8. 
In addition, the experimental values of the energies for 
fission isomers [103] are also shown by short horizontal 
lines. For 744Pu, a suggested isomeric energy range (1.6 
MeV to 2.6 MeV) from systematics was shown. From 
Fig. 8 and Table I one finds that for most of the nuclei 
investigated here, the triaxiality lowers the second barrier 
by about 0.5 ~ 1 MeV (about 10 ~ 20% of the barrier 
height) [46]. The calculation results with the triaxiality 
included agrees well with the empirical values for all ac- 
tinide nuclei shown in Fig. 8. In Ref. [46], we have found 
that for 24°Cm, the height of the second barrier without 
the triaxial deformation included is already smaller than 
the empirical value. From Table I it is seen that this is 
also the case for 7°°Cm and ?50:252Cf, For these four nu- 
clei, it is difficult to find the second saddle point when the 
triaxial deformation is included. Therefore the heights of 
the second barriers with the triaxial deformation included 
are not listed for them. The reason for these discrepan- 
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(Color online) Potential energy curves of even-even actinide nuclei around the ground states and the first fission 


barriers from MDC-RMF calculations. The axially symmetric results are displayed by dotted curves, while those from the 
triaxial calculations are shown by solid curves. The binding energy is normalized with respect to that of the ground state of 
each nucleus. The empirical values of fission barriers are taken from Ref. [102] and shown by red dots. 


cies may be related to a strong competition between the 
two or among more fission paths and the assumption on 
the barrier width made when the empirical value of the 
barrier height is evaluated, as discussed in Ref. [46]. 


D. Parameter dependence of the effect of 
triaxiality around the second barrier 


From the above detailed studies of PES’s of actinide 
nuclei in the (29, G22, 830) space, it is clear that the 
triaxiality plays an important role upon the second fis- 
sion barriers of actinide nuclei. We have studied the 
parameter dependence of the effect of triaxiality on the 
second barrier height and some results can be found in 
Refs. [46-48]. In Table II are listed the second barrier 
heights of ?34U calculated with different parameter sets, 
including meson-exchange ones NL3* [104], NL-Z2 [105], 
DD-ME2 [106], and point-coupling ones PC-PK1 [92] and 
DD-PC1 [107]. One finds that the second barrier height 
may differ by a few MeV with different effective inter- 
actions, but in all cases the barriers are lowered consid- 
erably by the nonaxial deformations. For this specific 
nucleus, the lowering effect with different effective inter- 


actions is about 0.5 MeV to 0.9 MeV which are clearly 
larger than the possible errors from the basis truncations 
discussed in Sec. III A. 


E. Further discussions 


A systematic study of fission barriers has been per- 
formed for even-even superheavy nuclei with Z = 112 ~ 
120 by using CDFTs and the outer fission barriers are 
found to be considerably affected by the triaxiality [49]. 
However, in a recent work within the MM model, it was 
found that the influence of the triaxiality on the second 
fission barriers is negligible [108]. This raises an open 
problem: What are the main reasons for the different ef- 
fects of the triaxiality on the second fission barriers pre- 
dicted by the MM model and CDFTs? It would be inter- 
esting to make more investigations with different models, 
especially, with the nonrelativistic DFTs. 

One of the problems concerning the PES calculated 
from self-consistent approaches is that there may ex- 
ist unexpected discontinuities on PES’s [60]. This is 
mainly due to the complexity of multi-dimensional PES’s. 
When a high-dimensional PES is projected onto a low- 
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FIG. 8. (Color online) Potential energy curves of even-even actinide nuclei around the fission isomers and the second fission 
barriers from MDC-RMF calculations. The reflection asymmetric shapes are allowed in the calculations. The axially symmetric 
results are displayed by dotted curves, while those from the triaxial calculations are shown by solid curves. The binding energy 
is normalized with respect to that of the ground state of each nucleus. The empirical values of fission barriers are taken from 
Ref. [102] and shown by red dots. Experimental values of energies of fission isomers are taken from Ref. [103] and shown by 
green horizontal lines. For 744Pu, only an energy range was given for the isomer. 


dimensional one, the minimizing procedure incorporated 
implicitly in the self-consistent approaches may cause 
a sudden change in the total energy at some points in 
the low-dimensional deformation space. This may re- 
sult in spurious saddle points. In our MDC-RMF cal- 
culations, we have tried some ways to avoid disconti- 
nuities on the PES’s and to exclude spurious saddle 
points [46]. To be completely free of discontinuities or 
spurious saddle points, one certainly should carry out 
multidimensionally-constrained calculations with higher- 
order multipole deformations included. In Ref. [109] the 
authors analyzed the origin of the discontinuities and pro- 
posed a numerical method to identify them. It will be 
useful to implement this method in the present MDC- 
RMF calculations; this will be one of our topics in the 
future. 


Since many symmetries are broken in mean field cal- 
culations, quantum numbers related to these symmetries 
are no longer conserved [2]. For example, the total an- 
gular momentum or nuclear spin J is not a good quan- 
tum number when (29 is not zero; similarly, the pro- 
jection of the total angular momentum on the symme- 


try axis K is also not conserved if a nucleus is triaxi- 
ally deformed [110]. When the reflection symmetry is 
broken, the parity is not a good quantum number. In 
our three-dimensional constraint calculations, all these 
quantum numbers are not good ones. One may intro- 
duce techniques of projection to restore the broken sym- 
metries [2]. In recent years, the importance of various 
projections on the PES’s and fission barriers have been 
investigated [39, 88, 111, 112]. 


IV. SUMMARY 


We developed multidimensionally-constrained rela- 
tivistic mean field (MDC-RMF) models. In these mod- 
els, the nuclear shape is assumed to be invariant under 
the reversion of x and y axes, i.e., the intrinsic sym- 
metry group is V4 and all shape degrees of freedom 
Pru with even H (B20, Boa, B30, B30, Bao, ven) are in- 
cluded self-consistently. The RMF functional can be 
one of the following four forms: the meson exchange or 
point-coupling nucleon interactions combined with the 
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TABLE I. Heights of the first and second fission barriers of 
some even-even actinide nuclei from MDC-RMF calculations 
compared with the empirical values (“Emp”) which are taken 
from Ref. [102]. The calculation results with and without the 
axial symmetry imposed are denoted by “AS” (axially sym- 
metric) and “TA” (triaxial), respectively. Note that around 
the first barrier, the reflection symmetry is assumed and for 
the second barrier, the octupole deformation is included. The 
height is in MeV. 


ZN A First barrier Second barrier 


AS TA Emp AS TA Emp 

20Th 90 140 230 5.03 3.96 6.10 6.80 6.37 6.50 
232Th 90 142 232 4.94 4.12 5.80 6.70 6.18 6.70 
2320 92 140 232 5.71 4.81 4.90 6.20 5.64 5.40 

23477 92 142 234 6.15 5.09 4.80 6.20 5.55 5.50 

2360 92 144 236 6.40 5.11 5.00 6.15 5.31 5.70 

2380 92 146 238 6.54 5.03 6.30 6.20 5.42 5.50 

2400 92 148 240 6.58 4.96 6.10 6.38 5.43 5.80 

238Pu 94 144 238 7.72 5.96 5.60 6.05 5.56 5.10 
240Pu 94 146 240 7.98 5.92 6.05 6.24 5.60 5.20 
242Pu 94 148 242 8.05 5.77 5.85 6.43 5.74 5.10 
244Pu 94 150 244 7.85 5.40 5.70 6.26 5.49 4.90 
246Pu 94 152 246 7.37 4.76 5.40 5.84 4.96 5.30 
242Cm 96 146 242 8.80 6.49 6.65 5.72 4.85 5.00 
440m 96 148 244 9.04 6.34 6.18 5.90 4.88 5.10 
460m 96 150 246 8.89 5.84 6.00 5.40 4.62 4.80 


Nucleus 


480m 96 152 248 8.43 5.35 5.80 4.10 — 4.80 
250Cm 96 154 250 7.77 4.79 5.40 2.60 — 4.40 
250Cf 98 152 250 8.87 5.70 5.60 2.40 — 3.80 


252Cf 98 154 252 8.41 5.26 5.30 1.20 — 3.50 


TABLE II. Height of the second barrier of 7°*U calculated 
from MDC-RMF models with different parameter sets. The 
results with and without nonaxial deformations included and 
their difference are denoted by Bas, Bra, and AB, respec- 
tively. 


Parameter set 


Bas (MeV) 


Bra (MeV) AB (MeV) 


NL3* 7.04 6.85 0.69 
NL-Z2 4.83 3.91 0.92 
PC-PK1 6.20 5.55 0.65 
DD-ME2 8.19 7.51 0.68 
DD-PC1 6.13 5.64 0.49 


nonlinear or density-dependent couplings. We solve the 
Dirac equation in an axially deformed harmonic oscil- 
lator (ADHO) basis. The convergence of the calculated 
results against the basis truncation, the oscillator length, 
and the basis deformation is studied and it is shown that 
reasonably large ADHO basis is able to provide a desired 
accuracy. 

Three-dimensional potential energy surface in the (620, 
b22, 839) plane is obtained for 74°Pu and potential energy 
curves of actinide nuclei around the first and second fis- 
sion barriers are studied systematically. It is found that 
the triaxiality is crucial in determining the height of both 
the first and the second fission barriers. Taking 734U as 
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an example, we have studied the parameter dependence 
of the effects of triaxiality on the second barrier height 
and found that the height of the second barrier may dif- 
fer by a few MeV with different effective interactions, but 
in all cases the barriers are lowered considerably by the 
nonaxial deformations. 

We conclude that it is important to include both the 
nonaxial and the spatial reflection asymmetric shapes si- 
multaneously for the study of potential energy surfaces 
and fission barriers of actinide nuclei and of those in un- 
known mass regions such as, e.g., superheavy nuclei. 
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Appendix A: Matrix elements of the Dirac 
Hamiltonian 


In the ADHO basis, the Dirac equation (5) is trans- 
formed into a matrix eigenvalue problem, 


Slam (r)a) fe’ +Y lajo - pla’)g” = ef? 
S (alo - pla’) f2’ + S(a|M_(r)la") 92" = eg? 


with M+ =+M +V 85. 
In the cylindrical coordinate system the kinetic energy 
operator o - p reads, 


= z e~’? (ð, — ip 1p) 
A ( e"? (8p + ip ðo) -ð; l 
The matrix elements of the kinetic energy operator o - p 
can be calculated as 


(alo : pla’) = (nz, np, Mi, Ms|a + pinz, np, My ms) 


ars 
= —1CQCa'bK,K! CE Sr, Aca 
/ 
+ ôm., 4 ôm, 1n ,n, (Baa — Mi Caa) 


2 
+ San. 10m! i Ôn (Baa +m}Caa’) ’ 
812 stg Varia 
(Al) 
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where the integrals read, 


Aga! = | di0n, 826, (A2) 
Bow = [ere ‘p0p ; (A3) 
Caw = J dpR™ RT, (A4) 


The matrix elements of the potentials U(r) = M4 (r) 
are calculated by expanding U(r) into a Fourier series, 


X UH (z, pe#?, 


u=- 


U(r) =U(z,p,9) = 5- (A5) 


where the components U) (z, p) are calculated by the 
inverse Fourier transformation, 


2m 
E | deU (z,p,p)e#, (A6) 
0 


where the integral over y is performed on a uniformly 
distributed mesh points. From the symmetry conditions 
one finds, 


UY) = yC”) = = (-1)4uU), (A7) 


Thus the Fourier series is abbreviated as 


U(r) = + (voe, p) +2 ` UC (z, p) cane) : 


n=1 
(A8) 
Finally the matrix element of U(r) reads, 
(a|U|a’) = = en ae ‘ pmp mi) 
= ôm, m, Cx Ca! == L [ör Di D®, 
+ 5 (0K'—K2n + ÔK-K',2n) D er J (A9) 
n=1 
where we have performed the integrals over y and 
Di = a a f pdpU (z, p)onebn, Rm RGA, 
(A10) 


are performed numerically by using the Gaussian quadra- 
ture. For axially symmetric potentials only the first term 
survives, while for triaxially deformed nucleus we must 
calculate additional terms with n Æ 0. 


Appendix B: Densities and their derivatives 


After solving the Dirac equation, we can calculate the 
densities from the wave functions. All the densities are 
linear combinations of the quantities pf (r, T) and p,(r,7) 
which are the contributions to the vector density from the 


13 


large and small components, respectively. 7 represents 
the isospin. For vector potential we have 


N 
eve Do we 


i=1 p=f.g 


5 (>: ve pepe Je T Oy. (B1) 


p=f,g ae’ \i=1 


pv(r) 


Note that the factors with y in ®, and ®, are in an 
exponential form, rearranging the terms, py(r) can be 
written as, 


pvr) = 5 : = (10) +25 eoon), 


n=1 
(B2) 
where the components 


N 
=>) >), (>: v pep? ) OR id esto 


p=f,g aa’ i=1 


XC2Ca! dn. Gn, Rie Bry | (B3) 
are calculated with the wave functions f% and gẹ. 

The scalar density ps(r) can be calculated similarly. 
The isovector densities are just the linear combinations 
of the corresponding densities of protons and neutrons. 

The derivatives of the vector density reads 


S (>. wpa Tia), (B4) 


p=f,g aa’ \i=1 
where the derivatives of the basis can be calculated as, 
VD ba) = Inet +2Jaar + Iža, (B5) 


with Iaa and Jaq defined later. First, from Eq. (9) we 
get 


V728,, = —2M [Ea — Va(z, p)| a, (B6) 

thus 
Toot = PIV? Dy 
= —2M |E. — Va(2, p)| msm, CC% 
mi l im 
Xn- Pn RR x —ellm—™)"_ (B7) 
a ERR eB 20 
Second, 
Jaw = VEl-VOq 


= Bray O Ca 
x fon, dnt (a, op Rma RA $ + mm) Fem Ri) 


m! 1 . 1 
+ RU Rp 0.bn,Oedn’ | x —eilm-m)e (B8) 
P Tig z 2T 
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Substituting them into Eq. (B5) and combining the same 
exponentials of y, we get 


1f. S (2n 
V’pv(r) = = (Wen +25 py” (z,p) cs(2ne) 
n=l 


[13] 
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